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Abstract: An iV-body code containing live stellar evolution through combination of the software 
packages NB0DY6 and STARS is presented. Operational details of the two codes are outlined and the 
changes that have been made to combine them discussed. We have computed the evolution of clusters 
of 10 4 stars using the combined code and we compare the results with those obtained using NB0DY6 
and the synthetic stellar evolution code SSE. We find that, providing the physics package within STARS 
is set up correctly to match the parameters of the models used to construct SSE, the results are very 
similar. This provides a good indication that the new code is working well. We also demonstrate 
how this physics can be changed simply in the new code with convective overshooting as an example. 
Similar changes in SSE would require considerable reworking of the model fits. We conclude by outlining 
proposed future development of the code to include more complete models of single stars and binary 
star systems. 
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1 Introduction 

The gravitational iV-body problem - the motion of a 
set of massive bodies under mutual self-gravity - is 
one of the oldest in computational astronomy . The 
first simulations were carried out by Tfolmbcrg ( 194ll ) 
who used special-purpose hardware. As computing 
power has increased, particul arly with the introduc - 
tion of the GRAPE hardware (|Makino and Taiiilll998l ') 
and computational algorithms have improved it has 
become possible to include greater numbers of parti- 
cles and more sophisticated physics in order to make 
more accurate models of star clusters, galaxies and the 
Universe itself. To study the evolution of galactic and 
globular open clusters we must integrate the motion of 
each particle separately, considering the force contri- 
bution from each of the other particles. Because of the 
high densities in the cores of stellar clusters particles 
can approach very close to one another and we must 
include special procedures to prevent such encounters 
introducing unacceptable errors. 

To make physically realistic models of stellar clus- 
ters we cannot avoid considering the evolution of the 
stars themselves. Mass loss from stars changes the 
cluster's potential and hence affects the motion of its 
stars. Furthermore the finite radius of a star means 
that if it approaches close to another they may inter- 
act. Tidal dissipation may cause a bound binary sys- 
tem to form, within which matter can be transferred 
from one star to another, and stars that come suffi- 
ciently close may merge. To predict the evolution of 
the mass and radius of a star with time it is necessary 



to use a stellar evolution code. Stellar evolution is 
another area of computational astronomy with a con- 
siderable heritage. A procedure for evaluating numer- 
ical models of st ars by use of an electr o nic co mputer 
was outlined by IHaselgrove and Hovld (|!956h . Over 
the years such simulations have grown in detail and 
accuracy as more physics has been added and more 
powerful computers have become available for the in- 
tegration of the equations. In the past however the 
computational cost of these calculations has prohib- 
ited including them in cluster simulations. 

The most popular approach to date for including 
stellar evolution in cluster models has been the use 
of synthetic stellar evolution codes. This involves an- 
alytic fits made to the results of a full code evalu- 
ated to simulate the evolu tion of a star . A pr ime ex- 
ample is the SSE code of iHurlev et al.1 (120021) which 
was d erived from the detailed models of iPols et al.l 
(1998). The combination of synthetic stellar evolution 
and iV-body dynamics has been used to good effect 
to model the interplay be tween stellar and dynamica l 
evolution in star c luster s llPortegies Zwart et al. I l200ll : 
iBaumgardt et al.1 [20031 : IHurlev et al.ll2005l ). However, 
a major drawback is that the fitting process is labori- 
ous and the result inflexible, in that it relies on choices 
made in generating the underlying set of detailed mod- 
els. If the input physics used for the detailed mod- 
els becomes out of date in any way then so does the 
synthetic package and it is non-trivial to generate a 
new set of analytic fits. It is also ill-equipped to deal 
with non-equilibrium cases arising from binary mass- 
transfer and stellar collisions. The pioneers of the syn- 



1 



2 



Publications of the Astronomical Society of Australia 



thetic stellar evolution approach were Eg gleton et al.l 
(1989) who preferred this to the main alternative at 
the time which was interpolation within a database of 
detailed models. It was decided that such a database 
would be cumbersome - especially if a large grid of 
models of various mass and metallicity were required 
- and create problems for data storage. Fortunately 
this is not a major issue anymore and the interpolation 
approach does have its merits. The ideal is to perform 
live stellar evolution in combination with JV-body dy- 
namics and to also include a module for dealing with 
any stellar collisions that a rise, for example t he smooth 
particle hydrodynamics of Sill s et all (|2003t ). 

Computing power has increased in recent times to 
the point where it is possible to combine a full stellar 
evolution code with a iV-body code. We present here 
such a code and compare the results to those obtained 
with synthetic stellar evolution. Section [5] contains de- 
tails of the two codes that we used and the models that 
we compare are described in Section[3] The results are 
presented in Section and discussed in Section [3] 

2 Description of the code 

We give a brief description of some of the relevant fea- 
tures of the two codes, STARS and NB0DY6, along with 
the modifications that have been necessary to combine 
them. 



the luminosity. A larger value of A allows the vari- 
ables to change more in a single timestep and hence 
larger timesteps are taken. Convective overshooting 
can be included by means of a modified Schwarzschild 
criterion controlled by a single param eter S ov ; details 
of the implementation are described in Schroder et al. 
119971 ). 

2.1.1 Zero-age models 

The models described in this paper were all calculated 
at solar composition, X = 0.7, Y = 0.28 and Z = 0.02. 
A library of zero-age main-sequence models of masses 
between O.IMq and 100 Mq has been produced by a 
multi-stage process. We took a stellar model of uni- 
form solar composition and inserted an artificial energy 
source to inflate it into a cold, low-density cloud until 
the core temperature was below 10 6 K, below which 
no nuclear reactions are modelled by the code. We 
then added and removed mass to produce a pre-main 
sequence of models of low density and temperature. 
We removed the artificial energy source and followed 
the contraction down to the point of minimum radius 
to produce a set of models with self-consistent central 
compositions and temperatures. By adding a small 
amount of mass to the envelope of one of these mod- 
els we can construct models of zero-age main-sequence 
stars with any mass. 



2.1 STARS 

STARS, the Camb ridge Stellar Evo lution code, was orig- 
inally written by lEggletonl ll 19711) and has b een exten- 
sively modified since (see iPols et all Il995l and refer- 
ences therein fo r a complete de s cripti on) . STARS utilises 
the method of iHenvev et al.l (|1964| ) . The equations 
of stellar structure are written as implicit finite dif- 
ference equations on a mesh and solved by numeri- 
cal inversion of the resulting matrix. The mesh has 
a fixed number of meshpoints that can move in both 
mass and radius, and hence is neither Eulerian nor La- 
grangian. The model is written in terms of the quan- 
tities log / (a measure of electr on degeneracy as de- 
scribed bv lEggleton et aT1ll973h . temperature, mass, 
radius, luminosity, the mass fractions of hydrogen, helium- 
4, carbon-12, oxygen-16 and neon-20 and a quantity Q 
which determines the position of the mesh. In a con- 
verged model the gradient of Q is equal at all points on 
the mesh. It has an ad-hoc functional form that causes 
more points to be placed where the temperature, mass, 
pressure and radius are varying most quickly. These 
are the regions, such as burning shells and ionisation 
zones in the atmosphere, that need to be well resolved. 
The system of equations is then solved by a relaxation 
method. The timestep in STARS is controlled by the 
user supplied parameter A. After a model has con- 
verged the next timestep is calculated as 



Sr i+ i — 5n 



A 



\sx jh \ 



(i) 



where 5Xjk is the change in variable j at meshpoint 
k. The sum is evaluated over all the variables except 



2.1.2 The helium flash 

We have developed a routine to pseudo-evolve through 
the helium flash in low-mass stars. In stars less mas- 
sive than about 2.3 Mq the core is degenerate at the 
time of helium ignition so that the increased tempera- 
ture owing to helium burning does not cause expansion 
and t hermal runaway occurs f Schwarzs child and Harml 
Il962f ). To evolve through the helium flash with STARS re- 
quires very small timesteps, which lead to numerical 
instability in calculation of the luminosity. To circum- 
vent these problems we construct approximate post- 
flash models with stable core helium burning. A star of 
mass M ~ 3 Mq that has evolved successfully through 
non-degenerate core helium ignition is taken and mat- 
ter removed from the envelope until the desired mass 
is reached. The hydrogen burning shell is then allowed 
to burn outwards with helium consumption disabled in 
order to obtain the correct core mass and the envelope 
compositions reset to their pre- flash values. Normal 
evolution is then resumed. Whilst not physically rig- 
orous this process provides models that can be used to 
simulate subsequent evolution. 



2.1.3 Mass loss 

We have written a simple procedure to determine the 
type of a star (main sequence, red giant, etc.) from the 
central abundances and hydrogen, helium and carbon 
burning luminosities of the star. Then we automati- 
cally choose an appropriate mass-loss law for the star 
according to 

• Main sequence - no mass loss, 
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Red Riant branch (RGB) - R eimers mass loss 
(|Kudritzki and Reimersll 1978ft 

with r) = 0.4, 

Core helium burning - Reimers mass loss with 
i? = 1.0, 

Asym ptotic giant branch (AGB) - Vassiliadi s and 
(|1993| ) mass loss 



AM 
~dT 



| 10 -11.4+0.0123P/d 3 ys MQyr -l 

L/c 



(-13.5 + 0.056P/days)kms~ 



(3) 



and 

• White dwarf (WD) - No mass loss. 

2.1.4 Timestep and mesh size 

To reduce model runtime and memory requirements 
we use a relatively low resolution (199 meshpoints per 
model). Runtime is further decreased by choosing a 
comparatively large value of the timestep parameter, 
A = 5, which ensures rapid progress throughout the 
star's life. These two choices have the added advantage 
of suppressing thermal pulses on the AGB, which are 
too computationally demanding and numerically dif- 
ficult to be included in these models (Stancliffe et al.l 
120041 ). 

2.1.5 The post-AGB 

Some extra procedures are necessary to complete the 
evolution of the stars from the late AGB to the white 
dwarf cooling track. This phase of evolution, known 
as the post-AGB, is even more numerically demand- 
ing than the AGB phase that proceeds it. To prevent 
numerical convergence issues arising from very high lu- 
minosities, once the stars have entered the superwind 
phase the rates of the triple-alpha and CNO reactions 
are capped at 10 -12 and 10 -11 reactions per baryon 
per second. Unmodified, STARS experiences a cyclical 
phase of evolution where fierce burning in the shell 
causes the very light envelope to be blown off. The 
shell extinguishes and then reignites as the envelope 
contracts back on to the star. This is both slow to con- 
verge, as many hundreds of cycles can occur even in the 
presence of strong mass loss and frequently causes code 
non-convergence. This behaviour is thought to be a 
numeric al artefact of th e STARS code though its cause is 
unclear (Stancliffe 2005). In order to surpress these cy- 
cles, once the hydrogen-burning shell reaches 0.02 Mq 
from the surface of the star the rate of energy gen- 
eration from nuclear burning is further reduced such 
that 



dE 
dt 



dE 



(4) 

v 0.02M o / dt physical 

where m cnv is the envelope mass (mass between the 
shell and surface of the star). The rate of consump- 
tion of material remains at the physical value, however, 



so the overall effect is to make the remaining nuclear 
burning produce less energy. 

As these procedures are both implemented only af- 
ter the superwind phase of mass loss has begun, which 
truncates the AGB, they have very little effect on the 
observable evolution. The most apparent consequence 
is that the luminosity on the post-AGB is reduced 
by about 0.5 dex; however, this phase is very short- 
liv ed and hence unlikely to be observed; no post-AGB 
Woffldaxs appear in the HR diagrams plotted Figures 1 to 
3. The principal effect of these two additional proce- 
dures is to render the predicted surface compositions 
of white dwarfs unreliable. These compositions are not 
trustworthy anyway because we do not model thermal 
pulses, which radically alter the compositions of AGB 
stars. We have found this procedure sufficient to evolve 
stars of initial masses up to 4Mq for a few Gyr, the 
typical lifetime of the clusters under consideration. 

2.1.6 Binaries 

The only element of binary stellar evolution imple- 
mented at present is the collision of stars. This is han- 
dled in a very simple manner. A zero- age star of the 
same mass as the combined mass of the two colliding 
stars is generated. Other interactions within binaries 
that form dynamically are ignored. 



2.2 NB0DY6 

NB0DY6 is a general-purpose full force summation N- 
body dynamics code. A general description of the de- 
velopment o f MB0DY6 and its sister codes can be found 
in a review ljAarsethlll9"99l ) and an exhaustive descrip- 
tion of its algorithmic b asis, construct ion and oper- 
ation in a recent book llAarsethl 120031 ). It uses the 
I Ahmad and Cohenl (|l973h neighbour scheme to reduce 
the cost of computation for large models, where the 
forces from particles lying within some local neigh- 
bour radius are updated more frequently than those 
from particles lying further away. An Hermite integra- 
tor is used with individual hierarchical timesteps. To 
remove the accuracy and performance problems associ- 
ated with integrating perturbed binarie s around many 
orbits Kustaanhcim o and Stiefel (|l965l ) regularisation 
is used. This involves a change of variables in four di- 
mensions, with one fictitious dimension added to make 
the transformation possible, as well as a time scal- 
ing. For perturbed hierarchical config urations (triples, 

quad ruples, etc.) chain regularisation ( Mik kola and Aarsethl 
1 19901 ) is used. By default MB0DY6 uses the s ynthetic 
stellar evolution code SSE (jHurlev et al.ll2000l ). 

To integrate STARS with NB0DY6 the procedures 
that involve calls to SSE have been re-written and a 
set of subroutines added to provide an interface be- 
tween the two codes. To interrogate the state of star 
Af at time t a routine has been written that extracts 
the required quantities from the live stellar evolution 
models at the time t (the evolution and A^-body codes 
have independent time steps). As a star evolves a 
model may occasionally fail to converge; in that case 
STARS restarts from the previous time with a reduced 
timestep. Hence it is possible for values at the cur- 
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rent timestep to change, because the restarted model 
may have different evolution. To avoid information 
provided to the A-body code being subsequently re- 
vised we ensure that the time for which the data are 
requested is between the previous and anteprevious 
timestep, i.e. that 

Wi-2 < t < TMi-1 

where tm% is the latest stellar evolution time for star 
M . If this is not the case then the model for star J\T 
is loaded into memory and evolved for the necessary 
number of timesteps. Once this has been done the 
radius, mass and luminosity at time t are calculated 
by linear interpolation between adjacent stellar mod- 
els and returned. A time for the next interrogation of 
the stellar model is also provided to NB0DY6. It is calcu- 
lated as the maximum time in which the radius should 
not change by more than 2% or the mass by more than 
1%. An arbitrary limit of ten STARS timesteps is also 
imposed, preventing NB0DY6 advancing past any inter- 
esting phases of evolution that start suddenly. 



3 Model setup 



To investigate the differences between models with syn- 
thetic and full stellar evolution five cluster models, 
each containing 10 4 stars, were evolved with our com- 
bined code. For the purposes of comparison we also 
evolved two clusters with the same initial conditions 
but with synthetic stellar evolution. Subsequent anal- 
ysis of the models produced caused us to run another 
model with full stellar evolution and convective over- 
shooting, with the parameter <5 v = 0.12 as in the con- 
struction of the SSE models. All stellar models started 
from the zero-age main sequence and had metallicity 
Z = 0.02. 

The stars are all taken to be single; there are no 
primordial binaries. The m asses of the st a rs are dis- 
tributed at random from a iKroupa et alJ (|1993l ) ini- 
tial mass function (IMF) , obtained from the generating 
function 



m(X) = 0.08 + 



7i A 72 + 73 
(1 - A) ' 58 



(5) 



where A is a number randomly chosen from the uni- 
form distribution in the range [0,1], 71 = 0.19, 72 = 
1.55, 73 = 0.05 and 74 = 0.6. The initial positions of 
the particles are distributed according to the Plummer 
distribution, 



p(r) = 



3M 



Anrl [1 + (r/r ) 2 ] 5 / 2 ' 



(6) 



where M is the total cluster mass and ro is a scal- 
ing factor related through integration to the half-mass 
radius rn by rn — 1.3 ro. Following the standard ap- 
proach for TV-body units and initial conditions we set 
M = 1 and ro = 1. The distance scaling of the sim- 
ulation is then determined by specifying the physical 
extent of the A-body length unit. We choose this to be 
2 pc for all these simulations which, along with the IMF 
and imposition of initial virial equilibrium, defines the 



scaling. More details can be found in I Aarse th (2003). 
A sta ndard tidal field based on Oort's constants (|Oortl 
1927) is applied which places the cluster at Sun's po- 
sition in the Galaxy. We take Oort's constants to be 
A = 14.4km s^kpc" 1 and B = -12.0km s" 1 kpc" 1 . 
The initial conditions for the five models were identi- 
cal except for the random number generator seed which 
was changed to give a different set of masses, positions 
and velocities for the stars. Hence the models can be 
considered to be a small, representative sample of pos- 
sible 10 4 star clusters with that set of initial conditions. 



4 Results and comparison 

We extracted various different properties from the data 
output by the code, chosen to measure different aspects 
of the stellar evolution and dynamics of the clusters in 
an attempt to determine how much the simulations 
differ. We have considered the quantities as a func- 
tion of the fraction of stars remaining instead of time 
because this is more characteristic of the dynamical 
state of evolution of the cluster and hence there is less 
scatter between the lines. In the case of the graph 
of the time against fraction of stars remaining this is 
counterintuitive but aids comparison with the other 
graphs in the paper. For most of the quantities we 
have also plotted the standard deviation of the values 
from the five standard models (full stellar evolution 
without convective overshooting). The points when 
the models have a certain fraction of stars remaining 
are not exactly coincident. Consequently we have in- 
terpolated the quantities to a standard set of points. 
This is a reasonable approach to take for quantities 
that are varying slowly, i.e. on a timescale longer than 
the interval between the snapshots of the cluster that 
the code produces. This standard deviation provides 
a measure of the spread of values owing to the varia- 
tion within the population of clusters of the type that 
we are considering. Furthermore we have calculated 
the mean of the five standard models at each of these 
points and subtracted it from the models using SSE 
and STARS with convective overshooting to get a set 
of residuals. By comparing these residuals with the 
standard deviations we measure whether the different 
stellar evolution prescriptions significantly change the 
models. 



4.1 HR diagrams 

The HR diagram of an actual cluster is relatively eas- 
ily observed - only photometry and the distance to the 
cluster are required for absolute magnitudes. Similarly 
the theoretical version, the luminosity-temperature di- 
agram, is easily generated as the temperature and lu- 
minosity are always calculated by a stellar evolution 
code. Comparison of the two requires bolometric cor- 
rection via fits to detailed models of stellar atmospheres. 
However, even given the uncertainties in this process, 
such diagrams provide a powerful and widely-used tool 
for comparing models and observations of clusters. Hence 
producing an accurate HR diagram is an important 
function of the stellar evolution section of a hybrid 
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code. Figure [T] contains HR diagrams for snapshots 
from one of the STARS models without convective over- 
shooting, Figure[2]HR diagrams for one of the SSE mod- 
els and Figure [3] diagrams for the STARS model com- 
puted with convective overshooting. 

4.2 Stellar type fractions and mass 

At the beginning of a model cluster's life its constituent 
stars are all on the zero-age main sequence. As time 
passes stars both evolve and are lost from the clus- 
ter. Both these processes affect the fraction of stars 
that are of a given type and hence these properties 
depend on both dynamics and stellar evolution. Fig- 
ure [J] contains diagrams showing the evolution of the 
various stellar type fractions. The same processes also 
affect the average stellar mass, the evolution of which 
is shown in Figure [5] The conversion of the fraction of 
stars remaining to time can be found from Figure [6] 

4.3 Time against number of stars 

Stars escape from a cluster by evaporation. Few-body 
interactions exchange energy between stars and an energy- 
gaining star can become unbound and escape from the 
cluster. This process is accelerated by the presence 
of the Galactic tidal field. Simultaneously the clus- 
ter loses mass because of evaporation, stellar winds 
and supernovae. This reduces the strength of the clus- 
ter potential and hence increases the evaporation rate. 
Thus the cluster mass and the number of stars remain- 
ing are tracers of dynamical processes happening on a 
smaller scale. A plot of the time taken for the fraction 
of stars remaining to fall to a given level is shown in 
Figure |S] and a plot of the fraction of the initial mass 
remaining in Figure [7] 

4.4 Core and half-mass radii 

The core is the most important part of a stellar cluster. 
In the core densities are highest, so the closest encoun- 
ters between stars take place. The crossing time of the 
core is smaller than that of the whole cluster, so catas- 
trophic collapse owing to the gravothermal instability 
takes place there first. The core radius is the stan- 
dard indicator of how large and how dense the core is, 
whilst the half-mass radius allows us to measure the 
behaviour of the outer parts of the cluster. A plot of 
the evolution of the core and half-mass radii with time 
is presented in Figure [5] 

5 Comparison and discussion 

The results show that the behaviour of the different 
codes is very similar, with some small differences which 
we discuss individually. 

5.1 Stellar measures 

Comparing Figures [T] and [5] several differences can be 
seen. The Hertzsprung gap in the STARS model is much 
more populated - this is particularly evident at 2 Gyr. 



The turnoff mass is lower in the STARS models because 
their main-sequence lifetimes are shorter. Finally the 
minimum white dwarf temperatures in the STARS mod- 
els are higher suggesting that they are cooling more 
slowly. This occurs even though the turnoff mass is 
lower, which implies that the white dwarfs formed ear- 
lier. 

The first two differences are caused by the presence 
of convective overshooting in the models from which 
SSE was constructed. Convective overshooting - the 
mixing of material outside the boundaries of convective 
regions in stars - is a possible explanation for the extra 
mixing over and above normal co nvective mixing that 
it tho ught to take place in stars ijShaviv and Salpeterl 
I1973T ). Amongst t he several effec t s it h as on evolu- 
tion, described in ISchroder et all (11997D . it extends 
the main-sequence lifetime and reduces the time spent 
crossing the Hertzsprung gap. The set of HR diagrams 
in Figure O for a cluster run using STARS with convec- 
tive overshooting, show that the first two differences 
largely disappear. Although there are still one or two 
more stars in the Hertzsprung gap the difference is not 
substantial. The differences in the white dwarf cool- 
ing track are caused by the over-simplified physics in- 
cluded in the version of SSE that we were using. This 
has been remedied in more recent versions. 

The overall trend in the evolution of the stellar 
type distributions, as shown in Figure [4] is that the 
fraction of main-sequence stars decreases as the clus- 
ter evolves whereas the fractions of white dwarfs and 
evolved stars increase. The decline in main-sequence 
star numbers is caused by stars evolving off the main 
sequence and the preferential evaporation of low-mass 
stars which are predominantly on the main sequence. 
The number of white dwarfs is larger than the num- 
ber of evolved stars other than at very early times 
because the evolved stars subsequently evolve further 
into white dwarfs. On the other hand as the cluster 
approaches dissolution the fraction of evolved stars in- 
creases much more rapidly, because these are now some 
of the heaviest stars in the cluster and hence the last 
to be lost. 

Again it can be seen from the plots of residuals 
that the inclusion of convective overshooting has a sig- 
nificant effect on cluster evolution. The fraction of 
main-sequence stars in the two models with convec- 
tive overshooting (the SSE models and the STARS model 
with convective overshooting) decline less quickly than 
in the STARS model without convective overshooting, 
particularly at the expense of evolved stars. For stars 
of intermediate mass overshooting enlarges the convec- 
tive core causing more hydrogen to be burnt and hence 
increasing the main-sequence lifetime. The SSE mod- 
els and the STARS model with convective overshooting 
are seen to agree within the intrinsic variability of the 
problem given by the error bars. 

5.2 Dynamical properties 

Comparing the results for SSE and STARS in Figure|H]we 
can see that for most of the models the total number 
of stars and mass of stars in the cluster are in good 
agreement between the two codes. The SSE models 
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Figure 1: HR diagrams for one of the STARS models with no convective overshooting. Different symbols 
represent different types of stars. The thick line of dots is the main sequence, and the sparser line of 
dots is the white dwarf cooling track. Open circles are Hertzsprung gap and first giant branch stars, and 
closed circles AGB stars. Filled squares are core helium burning stars. All binaries are assumed to be 
fully resolved (i.e. both stars are plotted separately). The HR diagrams are plotted at 500 Myr (top left), 
1000 Myr (top right), 2000 Myr (centre left), 3000Myr (centre right) and 4000Myr (bottom). The small 
numbers of stars where the model failed to converge - mostly stars with masses greater than 4 Mq, formed 
by collisions in perturbed binaries that have formed dynamically - have been removed from the diagrams. 
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Figure 2: As in Figure [T] for one of the SSE models. 
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Figure 3: As in Figured] for the STARS model with convective overshooting. 
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Figure 4: The top left panel shows the fraction of stars that are on the main sequence as a function of 
the fraction of stars remaining in the cluster. The solid lines are the STARS models without convective 
overshooting, the dashed lines the SSE models and the dotted line the STARS model with convective 
overshooting. The bottom left panel shows the deviations of the latter two values from the mean value 
of the five STARS models run without convective overshooting. The grey error bars show the standard 
deviation of these five models, the dashed lines the SSE models and the dotted line the STARS model with 
convective overshooting. The central panels show the same measurements for evolved stars (giants and 
core helium burning stars) and the right panels those for white dwarfs. 
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Figure 5: The left panel shows the average stellar mass in the models plotted against the fraction of 
stars remaining in the cluster. The right panel shows the deviations of the values from the mean of 
the five STARS models calculated without convective overshooting. The solid lines are the STARS models 
without convective overshooting, the dashed lines the SSE models and the dotted line the STARS model 
with convective overshooting. The error bars represent the standard deviation of the five STARS models 
without convective overshooting. 




N/N„ N/N 



Figure 6: As in Figure [5j but for the physical time as a function of the fraction of stars remaining in the 
cluster. 
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Figure 8: The evolution of core radius (lower, grey line) and half-mass radius (upper, black line) with 
fraction of stars remaining. The values of the two radii have been smoothed by averaging across three 
consecutive values. The five solid lines are the STARS models without convective overshooting, the dashed 
line one of the SSE models and the dotted line the STARS model with convective overshooting. The right 
panel shows just the first part of the evolution to facilitate observation of core collapse. 
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have fewer stars towards the end of the evolutionary 
sequence (after about half the stars have been lost from 
the simulations). To explain this difference it is useful 
to consider first the average stellar mass. 

We can see in Figure [5] that the average stellar 
mass behaves similarly across all the models. Initially 
it declines slightly owing to mass loss from the most 
massive stars and then increases following the prefer- 
ential evaporation of low-mass stars. From the plot of 
the residuals we see that until about half the stars in 
the cluster have been lost the SSE models have higher 
average masses. This is because the inclusion of con- 
vective overshooting increases the main-sequence life- 
times of intermediate-mass stars. This in turn in- 
creases the turnoff mass and hence the mean stellar 
mass. However as the cluster evolves the stars at 
the turnoff have smaller and smaller convective cores 
and so the difference in lifetimes is less pronounced. 
The residuals of the average mass for the STARS model 
with convective overshooting shows a similar trend but 
displaced downward slightly, as the mean mass drops 
more sharply at the beginning of the simulation. This 
is presumably because there are, by chance, a larger 
number of the highest mass stars in the simulation. 
These effects are also clear in the plot of the total mass 
remaining in the cluster (Figure 0. 

The effect of this enhanced number of more mas- 
sive stars appears to be to increase the rate of ejection 
of stars from the cluster. This explains the number 
of stars being lower in the SSE models from the point 
where half the stars have evaporated onwards. How- 
ever, one should be aware of reading too much into 
the results; few of the differences are significant be- 
yond two standard deviations. 

From the plot of core and half- mass radii (Figure[5)l 
it is hard to glean much interesting information. Look- 
ing at the half-mass radius we can see that the clus- 
ter expands initially, owing to mass loss from evolving 
stars, then remains at roughly constant size for most 
of its lifetime and shrinks towards the end as it dis- 
solves. The core radius declines sharply at the start 
of the simulation, undergoes core collapse just prior 
to 200 Myr and behaves erratically thereafter. No sig- 
nificant deviation between the different sets of stellar 
models can be identified. 



6 Conclusions 

We would expect, a priori, that the introduction of 
full stellar evolution would have a limited effect on 
cluster simulations containing only single stars because 
the synthetic stellar evolution fits are generally found 
to be good for single stars. These results show that 
this is indeed the case. Once convective overshooting 
has been included in the models they fit the results 
obtained with SSE to within the intrinsic variability of 
the calculations. 

One of the benefits of a live stellar evolution code 
is the increased flexibility that it brings. In order to 
change the value of the convective overshoot parame- 
ter, the mixing length parameter for convective mix- 
ing, the initial elemental abundances, or any other as- 



pect of the physics package, inside STARS we merely a 
change one or two parameters. This contrasts with a 
synthetic code where a whole new grid of models must 
be computed and fits made to them. Hence a live code, 
whilst not offering much improvement in accuracy for 
single stars, does provide additional flexibility. The 
new code allows us to produce, for instance, stellar 
isochrones and luminosity functions that incorporate 
the effect of dynamics, whilst also allowing us to vary 
aspects of the stellar physics such as the mixing and 
initial chemical compositions. 

To make full use of detailed stellar evolution mod- 
els we need to extend the code described here. The 
problems with numerical non-convergence on the late 
AGB/post-AGB need to be solved for stars of mass 
greater than 4Mq so that we can simulate the whole 
range of masses of single stars. The next step will then 
be to add the effects of binary interactions. Interac- 
tions are very significant for the dynamical evolution 
of clusters and binaries are responsible for the forma- 
tion of many different types of stellar exotica, includ- 
ing blue stragglers, X-ray binaries, millisecond pulsars, 
etc. For a full description of cluster evolution, this im- 
portant (and often difficult) area of evolution needs to 
be modelled properly. Such processes as Roche Lobe 
overflow and common envelope evolution must be in- 
cluded. Work on developing the code further to this 
end is underway. 
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